This tutorial is part of multiphenicsx.

Tutorial 01. Block Poisson problem.

In this tutorial we first solve the problem

\begin{cases} -u'' = f & \text{in }\Omega = (0, 1),\\ u = 0 & \text{on }\Gamma = \{0, 1\} \end{cases}

using standard (non block) dolfinx code.

Then we use block support of dolfinx to solve the system

\begin{cases} - w_1'' - 2 w_2'' = 3 f & \text{in }\Omega,\\ - 3 w_1'' - 4 w_2'' = 7 f & \text{in }\Omega \end{cases}

subject to

\begin{cases} w_1 = 0 & \text{on }\Gamma,\\ w_2 = 0 & \text{on }\Gamma \end{cases}

By construction the solution of the system is $$(w_1, w_2) = (u, u)$$

We then compare the obtained solutions. This tutorial serves as a reminder on how to solve a linear problem in dolfinx, and introduces multiphenicsx.fem.petsc.BlockVecSubVectorWrapper.

import typing

import dolfinx.fem
import dolfinx.fem.petsc
import dolfinx.mesh
import mpi4py.MPI
import numpy as np
import numpy.typing
import petsc4py.PETSc
import ufl
import viskex

import multiphenicsx.fem
import multiphenicsx.fem.petsc

Mesh generation

mesh = dolfinx.mesh.create_unit_interval(mpi4py.MPI.COMM_WORLD, 32)

def left(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[np.bool_]:
    """Identify the left boundary of the domain."""
    return abs(x[0] - 0.) < np.finfo(float).eps  # type: ignore[no-any-return]


def right(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[np.bool_]:
    """Identify the right boundary of the domain."""
    return abs(x[0] - 1.) < np.finfo(float).eps  # type: ignore[no-any-return]


left_facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim - 1, left)
right_facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim - 1, right)
boundary_facets = np.hstack((left_facets, right_facets))

x0 = ufl.SpatialCoordinate(mesh)[0]

viskex.dolfinx.plot_mesh(mesh)

Standard case (solve for $u$)

def run_standard() -> dolfinx.fem.Function:
    """Run the standard case: solve for the unknown u."""
    # Define a function space
    V = dolfinx.fem.functionspace(mesh, ("Lagrange", 2))
    u = ufl.TrialFunction(V)
    v = ufl.TestFunction(V)

    # Define problems forms
    a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
    f = ufl.inner(100 * ufl.sin(20 * x0), v) * ufl.dx

    # Define boundary conditions
    zero = petsc4py.PETSc.ScalarType(0)
    bdofs_V = dolfinx.fem.locate_dofs_topological(V, mesh.topology.dim - 1, boundary_facets)
    bc = dolfinx.fem.dirichletbc(zero, bdofs_V, V)

    # Solve the linear system
    u = dolfinx.fem.Function(V)
    problem = dolfinx.fem.petsc.LinearProblem(
        a, f, bcs=[bc], u=u,
        petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})
    problem.solve()
    u.vector.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)

    # Return the solution
    return u  # type: ignore[no-any-return]

u = run_standard()

viskex.dolfinx.plot_scalar_field(u, "u")

Block case (solve for ($w_1$, $w_2$))

def run_block() -> typing.Tuple[dolfinx.fem.Function, dolfinx.fem.Function]:
    """Run the block case: solve for the unknowns (w_1, w_2)."""
    # Define a block function space
    V1 = dolfinx.fem.functionspace(mesh, ("Lagrange", 2))
    V2 = dolfinx.fem.functionspace(mesh, ("Lagrange", 2))
    (u1, u2) = (ufl.TrialFunction(V1), ufl.TrialFunction(V2))
    (v1, v2) = (ufl.TestFunction(V1), ufl.TestFunction(V2))

    # Define problem block forms
    aa = [[1 * ufl.inner(ufl.grad(u1), ufl.grad(v1)) * ufl.dx, 2 * ufl.inner(ufl.grad(u2), ufl.grad(v1)) * ufl.dx],
          [3 * ufl.inner(ufl.grad(u1), ufl.grad(v2)) * ufl.dx, 4 * ufl.inner(ufl.grad(u2), ufl.grad(v2)) * ufl.dx]]
    ff = [ufl.inner(300 * ufl.sin(20 * x0), v1) * ufl.dx,
          ufl.inner(700 * ufl.sin(20 * x0), v2) * ufl.dx]
    aa_cpp = dolfinx.fem.form(aa)
    ff_cpp = dolfinx.fem.form(ff)

    # Define block boundary conditions
    zero = petsc4py.PETSc.ScalarType(0)
    bdofs_V1 = dolfinx.fem.locate_dofs_topological(V1, mesh.topology.dim - 1, boundary_facets)
    bdofs_V2 = dolfinx.fem.locate_dofs_topological(V2, mesh.topology.dim - 1, boundary_facets)
    bc1 = dolfinx.fem.dirichletbc(zero, bdofs_V1, V1)
    bc2 = dolfinx.fem.dirichletbc(zero, bdofs_V2, V2)
    bcs = [bc1, bc2]

    # Assemble the block linear system
    AA = dolfinx.fem.petsc.assemble_matrix_block(aa_cpp, bcs)
    AA.assemble()
    FF = dolfinx.fem.petsc.assemble_vector_block(ff_cpp, aa_cpp, bcs)

    # Solve the block linear system
    uu = dolfinx.fem.petsc.create_vector_block(ff_cpp)
    ksp = petsc4py.PETSc.KSP()
    ksp.create(mesh.comm)
    ksp.setOperators(AA)
    ksp.setType("preonly")
    ksp.getPC().setType("lu")
    ksp.getPC().setFactorSolverType("mumps")
    ksp.setFromOptions()
    ksp.solve(FF, uu)
    uu.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)
    ksp.destroy()

    # Split the block solution in components
    u1u2 = (dolfinx.fem.Function(V1), dolfinx.fem.Function(V2))
    with multiphenicsx.fem.petsc.BlockVecSubVectorWrapper(uu, [V1.dofmap, V2.dofmap]) as uu_wrapper:
        for u_wrapper_local, component in zip(uu_wrapper, u1u2):
            with component.vector.localForm() as component_local:
                component_local[:] = u_wrapper_local

    # Return the block solution components
    return u1u2

u1, u2 = run_block()

viskex.dolfinx.plot_scalar_field(u1, "u1")

viskex.dolfinx.plot_scalar_field(u2, "u2")

Error computation between the different cases

def compute_errors(
    u1: dolfinx.fem.Function, u2: dolfinx.fem.Function, uu1: dolfinx.fem.Function, uu2: dolfinx.fem.Function
) -> None:
    """Compute errors between standard and standard block cases."""
    u_1_norm = np.sqrt(mesh.comm.allreduce(
        dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(ufl.grad(u1), ufl.grad(u1)) * ufl.dx)),
        op=mpi4py.MPI.SUM))
    u_2_norm = np.sqrt(mesh.comm.allreduce(
        dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(ufl.grad(u2), ufl.grad(u2)) * ufl.dx)),
        op=mpi4py.MPI.SUM))
    err_1_norm = np.sqrt(mesh.comm.allreduce(
        dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(ufl.grad(u1 - uu1), ufl.grad(u1 - uu1)) * ufl.dx)),
        op=mpi4py.MPI.SUM))
    err_2_norm = np.sqrt(mesh.comm.allreduce(
        dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(ufl.grad(u2 - uu2), ufl.grad(u2 - uu2)) * ufl.dx)),
        op=mpi4py.MPI.SUM))
    print("  Relative error for first component is equal to", err_1_norm / u_1_norm)
    print("  Relative error for second component is equal to", err_2_norm / u_2_norm)
    assert np.isclose(err_1_norm / u_1_norm, 0., atol=1.e-10)
    assert np.isclose(err_2_norm / u_2_norm, 0., atol=1.e-10)

print("Computing errors between standard and block")
compute_errors(u, u, u1, u2)

Computing errors between standard and block

  Relative error for first component is equal to 4.923782490298636e-15
  Relative error for second component is equal to 2.7591822783167812e-14